home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / bit / src / ulib / interpol.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  4KB  |  189 lines
  1. /***********************************************************************
  2.  * $Id: interpol.c,v 0.80 1994/02/24 09:48:11 zhao Exp $
  3.  *
  4.  *.  Copyright(c) 1993,1994 by T.C. Zhao
  5.  *   All rights reserved.
  6.  *.
  7.  *
  8.  * Interpolate a tabulated function onto a uniform grid using Nth order
  9.  * Lagrangian polynomial. deg specifies the degree of interpolating
  10.  * polynomial but in practice, deg > 3 is a bad idea.
  11.  *
  12.  * If we can assume input is uniform grided, Newton's forward difference
  13.  * polynomial method can be used, which is faster.
  14.  *
  15.  ***********************************************************************/
  16. #if !defined(lint) && defined(F_ID)
  17. char *id_intpl = "$Id: interpol.c,v 0.80 1994/02/24 09:48:11 zhao Exp $";
  18. #endif
  19.  
  20. #include "ulib.h"
  21.  
  22. int
  23. lp_interpol(float *xin, float *yin, int nin, float grid,
  24.         float *xout, float *yout, int *nout, int deg)
  25. {
  26.  
  27.     int i, j, k, l, n, jo;
  28.     int ih, im, idm;
  29.     float term;
  30.  
  31.     /* sanity checks */
  32.  
  33.     if (grid < 0.0 || nin < (deg + 1))
  34.       {
  35.       M_err("Interpol", "Bad arg grid=%g nin=%d", grid, nin);
  36.       return -1;
  37.       }
  38.  
  39.     /* points in output */
  40.     n = ((xin[nin - 1] - xin[0]) / grid + 0.01) + 1;
  41.     if (n > *nout)
  42.       {
  43.       M_warn("Interpol", "output truncated");
  44.       n = *nout;
  45.       }
  46.  
  47.     xout[0] = xin[0];
  48.     yout[0] = yin[0];
  49.  
  50.     jo = 0;
  51.     for (i = 1; i < n; i++)
  52.       {
  53.       xout[i] = xin[0] + i * grid;
  54.  
  55.       /* center current point using binary search */
  56.       j = jo;
  57.       ih = nin;
  58.       while ((ih - j) > 1)
  59.         {
  60.         im = (ih + j) / 2;
  61.         if (xout[i] > xin[im])
  62.             j = im;
  63.         else
  64.             ih = im;
  65.         }
  66.  
  67.       jo = j;
  68.  
  69.       j -= deg / 2;
  70.  
  71.       /* boundary checks */
  72.       if (j < 0)
  73.           j = 0;
  74.  
  75.       if (j > (nin - deg - 1))
  76.           j = nin - deg - 1;
  77.  
  78.       /* do it now */
  79.       yout[i] = 0.0;
  80.       idm = j + deg;
  81.  
  82.       for (l = j; l <= idm; l++)
  83.         {
  84.         term = yin[l];
  85.         for (k = j; k <= idm; k++)
  86.           {
  87.               if (l != k)
  88.               term *= (xout[i] - xin[k]) / (xin[l] - xin[k]);
  89.           }
  90.         yout[i] += term;
  91.         }
  92.       }
  93.     *nout = n;
  94.     return 0;
  95. }
  96.  
  97. /*********************************************************************
  98.  * A single point
  99.  ********************************************************************/
  100.  
  101. float
  102. val_lp_interpol(float *x, float *y, int n, float xval, int deg)
  103. {
  104.     float val, term;
  105.     int j, l, k;
  106.     int ih, im, idm;
  107.  
  108.     if (xval < x[0] || xval > x[n - 1])
  109.       {
  110.       M_err("ValInterpol", "Failed to bracket %g", x);
  111.       return 0.0;
  112.       }
  113.  
  114.     /* hunt to center */
  115.     j = 0;
  116.     ih = n;
  117.     while ((ih - j) > 1)
  118.       {
  119.       im = (ih + j) / 2;
  120.       if (xval > x[im])
  121.           j = im;
  122.       else
  123.           ih = im;
  124.       }
  125.  
  126.  
  127.     j -= deg / 2;
  128.  
  129.     /* boundary checks */
  130.     if (j < 0)
  131.     j = 0;
  132.  
  133.     if (j > (n - deg - 1))
  134.     j = n - deg - 1;
  135.  
  136.     /* do it now */
  137.     val = 0.0;
  138.     idm = j + deg;
  139.  
  140.     for (l = j; l <= idm; l++)
  141.       {
  142.       term = y[l];
  143.       for (k = j; k <= idm; k++)
  144.         {
  145.         if (l != k)
  146.             term *= (xval - x[k]) / (x[l] - x[k]);
  147.         }
  148.       val += term;
  149.       }
  150.     return val;
  151. }
  152.  
  153. #ifdef TEST
  154. #include <stdio.h>
  155.  
  156. main()
  157. {
  158.     float grid = 0.2;
  159.     float x[100], y[100];
  160.     float inx[30], iny[30];
  161.     float *xp = inx, *yp = iny;
  162.     int i = 0, nout;
  163.  
  164.  
  165.     while (scanf(" %f %f ", xp, yp) == 2)
  166.       {
  167.       xp++;
  168.       yp++;
  169.       i++;
  170.       }
  171.  
  172.     nout = 100;
  173.  
  174.     fprintf(stderr, "NIN=%d\n", i);
  175. #if 0
  176.     lp_interpol(inx, iny, i, grid, x, y, &nout, 3);
  177.     xp = x;
  178.     yp = y;
  179.     printf("N=%d\n", nout);
  180.     for (i = 0; i < nout; i++, xp++, yp++)
  181.     printf("%8.3f %8.3f\n", *xp, *yp);
  182. #else
  183.     printf("X=%f Y=%f\n",
  184.        3.32, val_lp_interpol(inx, iny, i, 3.32, 3));
  185. #endif
  186. }
  187.  
  188. #endif
  189.